GARCH in R

install.packages("PerformanceAnalytics", type = "source")
## Installing package into '/opt/homebrew/lib/R/4.4/site-library'
## (as 'lib' is unspecified)
library(xts)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(zoo)
library(PerformanceAnalytics)
## 
## Attaching package: 'PerformanceAnalytics'
## The following object is masked from 'package:graphics':
## 
##     legend
getwd()
## [1] "/Users/michal/Documents/Code/metals"
# setwd("C:/Users/user/OneDrive/Dokumenty/R/classes2024")
list.files()
##  [1] "__pycache__"              "bloomberg_data"          
##  [3] "combined_data.csv"        "data"                    
##  [5] "data_loading.ipynb"       "data_vis.ipynb"          
##  [7] "fetch_dailymetalprice.py" "GARCH_S&P500.ipynb"      
##  [9] "GARCH.html"               "GARCH.ipynb"             
## [11] "GARCH.Rmd"                "README.md"               
## [13] "sp500prCl.csv"            "utils"
setwd("/Users/michal/Documents/Code/metals")
sp500prices <- read.csv("sp500prCl.csv", dec = ".", sep = ";")
# sp500prices
# Plot daily S&P 500 prices
plot(sp500prices$Price)

class(sp500prices)
## [1] "data.frame"
stockdata <- sp500prices
date <- as.POSIXct(strptime(stockdata$Data, format = "%Y-%m-%d")) # "%Y-%m-%d"  "%d.%m.%Y"
# date
sp500prices1 <- xts(stockdata[2:ncol(stockdata)], order.by = date) #
rm(stockdata)

head(sp500prices1$Price)
##             Price
## 1989-01-03 275.31
## 1989-01-04 279.43
## 1989-01-05 280.01
## 1989-01-06 280.67
## 1989-01-09 280.98
## 1989-01-10 280.38
# Compute daily returns
# sp500ret <- diff(log(sp500prices1))
sp500ret <- CalculateReturns(sp500prices1) # watch the dog! first ret is zero!
sp500ret <- sp500ret[-1, ]
# Check the class of sp500ret
class(sp500ret)
## [1] "xts" "zoo"
# Plot daily returns
plot(sp500ret)

# Compute the daily standard deviation for the complete sample
sd(sp500ret)
## [1] 0.0109947
# Compute the annualized volatility for the complete sample
sqrt(252) * sd(sp500ret)
## [1] 0.1745354
# Compute the annualized standard deviation for the year 2009
sqrt(252) * sd(sp500ret["2009"])
## [1] 0.2728327
# Compute the annualized standard deviation for the year 2017
sqrt(252) * sd(sp500ret["2017"])
## [1] 0.06685658
# Showing two plots on the same figure
par(mfrow = c(2, 1))

# Compute the rolling 1 month estimate of annualized volatility
chart.RollingPerformance(
     R = sp500ret["2000::2017"], width = 22,
     FUN = "sd.annualized", scale = 252, main = "One month rolling volatility"
)

# Compute the rolling 3 months estimate of annualized volatility
chart.RollingPerformance(
     R = sp500ret["2000::2017"], width = 66,
     FUN = "sd.annualized", scale = 252, main = "Three months rolling volatility"
)

## Compute absolute value of the prediction errors

# Compute the mean daily return
m <- mean(sp500ret)

# Define the series of prediction errors
e <- sp500ret - m

# Plot the absolute value of the prediction errors
par(mfrow = c(2, 1), mar = c(3, 2, 2, 2))
plot(abs(e))

# Plot the acf of the absolute prediction errors
acf(abs(e), lag.max = 20)

# GARCH:

# Define parameters and numerics
alpha <- 0.1
beta <- 0.8
omega <- var(sp500ret) * (1 - alpha - beta)
e <- sp500ret - mean(sp500ret)
e2 <- e^2
nobs <- length(sp500ret)
predvar <- rep(NA, nobs)

# Compute the predicted variances
predvar[1] <- var(sp500ret)
for (t in 2:nobs) {
     predvar[t] <- omega + alpha * e2[t - 1] + beta * predvar[t - 1]
}

# Create annualized predicted volatility
ann_predvol <- xts(sqrt(252) * sqrt(predvar), order.by = time(sp500ret))

# Plot the annual predicted volatility in 2008 and 2009
plot(ann_predvol["2008::2009"], main = "Ann. S&P 500 vol in 2008-2009")

# Time for rugarch package

# install.packages("rugarch")
library(rugarch)
## Loading required package: parallel
## 
## Attaching package: 'rugarch'
## The following object is masked from 'package:stats':
## 
##     sigma
# Specify a standard GARCH model with constant mean
garchspec <- ugarchspec(
     mean.model = list(armaOrder = c(0, 0)),
     variance.model = list(model = "sGARCH"),
     distribution.model = "norm"
)

# Estimate the model
garchfit <- ugarchfit(data = sp500ret, spec = garchspec)

# See the resutls:
garchfit
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(1,1)
## Mean Model   : ARFIMA(0,0,0)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.000572    0.000087  6.59621 0.000000
## omega   0.000001    0.000002  0.67008 0.502805
## alpha1  0.078259    0.027893  2.80567 0.005021
## beta1   0.910553    0.029458 30.91016 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.000572    0.000853 0.670332  0.50265
## omega   0.000001    0.000058 0.021204  0.98308
## alpha1  0.078259    0.901416 0.086818  0.93082
## beta1   0.910553    0.946490 0.962031  0.33603
## 
## LogLikelihood : 24073.8 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.5881
## Bayes        -6.5844
## Shibata      -6.5881
## Hannan-Quinn -6.5869
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                      1.311 0.25213
## Lag[2*(p+q)+(p+q)-1][2]     1.582 0.34264
## Lag[4*(p+q)+(p+q)-1][5]     7.253 0.04502
## d.o.f=0
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                    0.07822  0.7797
## Lag[2*(p+q)+(p+q)-1][5]   5.15209  0.1414
## Lag[4*(p+q)+(p+q)-1][9]   6.34246  0.2607
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]    0.1801 0.500 2.000  0.6713
## ARCH Lag[5]    0.2836 1.440 1.667  0.9446
## ARCH Lag[7]    0.8098 2.315 1.543  0.9424
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  346.8562
## Individual Statistics:               
## mu      0.04305
## omega  30.47592
## alpha1  0.20103
## beta1   0.29761
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.07 1.24 1.6
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value      prob sig
## Sign Bias           2.6152 8.935e-03 ***
## Negative Sign Bias  0.3192 7.496e-01    
## Positive Sign Bias  2.3649 1.806e-02  **
## Joint Effect       33.3642 2.698e-07 ***
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     223.6    6.447e-37
## 2    30     232.5    1.196e-33
## 3    40     270.2    2.282e-36
## 4    50     293.3    1.537e-36
## 
## 
## Elapsed time : 0.2587059
# Use the method sigma to retrieve the estimated volatilities
garchvol <- sigma(garchfit)
garch_mean <- fitted(garchfit)

# Plot the volatility for 2017
plot(garchvol["2017"])

# Forecasts:

# Compute unconditional volatility
sqrt(uncvariance(garchfit))
## [1] 0.01051693
# Print last 10 ones in garchvol
tail(garchvol, 10)
##                   [,1]
## 2017-12-15 0.005071375
## 2017-12-18 0.005493771
## 2017-12-19 0.005524124
## 2017-12-20 0.005491357
## 2017-12-21 0.005371077
## 2017-12-22 0.005259472
## 2017-12-26 0.005148613
## 2017-12-27 0.005057916
## 2017-12-28 0.004953331
## 2017-12-29 0.004868582
# Forecast volatility 5 days ahead and add
garchforecast <- ugarchforecast(fitORspec = garchfit, n.ahead = 5)

# Extract the predicted volatilities and print them
print(sigma(garchforecast))
##      2017-12-29
## T+1 0.005041038
## T+2 0.005134710
## T+3 0.005225683
## T+4 0.005314106
## T+5 0.005400117

Volatility targeting

# When you target a portfolio with 5% annualized volatility, and the annualized volatility of the risky asset is, then you should invest in the risky asset.

# Compute the annualized volatility
annualvol <- sqrt(252) * sigma(garchfit)

# Compute the 5% vol target weights
vt_weights <- 0.05 / annualvol

# Compare the annualized volatility to the portfolio weights in a plot
plot(merge(annualvol, vt_weights), multi.panel = TRUE)

# Application of the GARCH models to the data:

# generate ret series"
# list(mu = 0, ar1 = 0, ma1 = 0, omega = 6*10^(-7), alpha1 = 0.07, beta1 = 0.9, skew = 0.9, shape = 5)
ret <- sp500ret
# Plot the return series

plot(ret)

# Specify the garch model to be used
garchspec <- ugarchspec(
     mean.model = list(armaOrder = c(0, 0)),
     variance.model = list(model = "sGARCH"),
     distribution.model = "sstd"
)

# Estimate the model
garchfit <- ugarchfit(data = ret, spec = garchspec)
garchvol <- sigma(garchfit)
garch_mean <- fitted(garchfit)
# Inspect the coefficients
coef(garchfit)
##           mu        omega       alpha1        beta1         skew        shape 
## 5.648052e-04 6.398380e-07 7.521339e-02 9.218343e-01 9.435684e-01 6.315897e+00
# See the whole output
garchfit
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(1,1)
## Mean Model   : ARFIMA(0,0,0)
## Distribution : sstd 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.000565    0.000087   6.4565 0.000000
## omega   0.000001    0.000000   2.7701 0.005604
## alpha1  0.075213    0.004718  15.9414 0.000000
## beta1   0.921834    0.004427 208.2420 0.000000
## skew    0.943568    0.014671  64.3136 0.000000
## shape   6.315897    0.454050  13.9101 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.000565    0.000081   6.9867 0.000000
## omega   0.000001    0.000001   1.0855 0.277682
## alpha1  0.075213    0.020742   3.6262 0.000288
## beta1   0.921834    0.018987  48.5506 0.000000
## skew    0.943568    0.014018  67.3091 0.000000
## shape   6.315897    0.627827  10.0599 0.000000
## 
## LogLikelihood : 24279.25 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.6438
## Bayes        -6.6382
## Shibata      -6.6438
## Hannan-Quinn -6.6419
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                      1.197 0.27384
## Lag[2*(p+q)+(p+q)-1][2]     1.498 0.36140
## Lag[4*(p+q)+(p+q)-1][5]     7.382 0.04186
## d.o.f=0
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                   0.002836  0.9575
## Lag[2*(p+q)+(p+q)-1][5]  4.322017  0.2165
## Lag[4*(p+q)+(p+q)-1][9]  5.564033  0.3514
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]    0.2192 0.500 2.000  0.6396
## ARCH Lag[5]    0.2926 1.440 1.667  0.9422
## ARCH Lag[7]    1.0704 2.315 1.543  0.9018
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  748.9782
## Individual Statistics:                
## mu       0.07006
## omega  156.38935
## alpha1   0.37783
## beta1    0.56251
## skew     0.41544
## shape    0.71577
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.49 1.68 2.12
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value      prob sig
## Sign Bias           2.7461 6.045e-03 ***
## Negative Sign Bias  0.1844 8.537e-01    
## Positive Sign Bias  2.5478 1.086e-02  **
## Joint Effect       33.9210 2.059e-07 ***
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     53.96    3.349e-05
## 2    30     61.00    4.632e-04
## 3    40     71.14    1.263e-03
## 4    50     83.73    1.463e-03
## 
## 
## Elapsed time : 0.38465
# Compute the standardized returns
stdret <- residuals(garchfit, standardize = TRUE)

# Compute the standardized returns using fitted() and sigma()
stdret <- (ret - fitted(garchfit)) / sigma(garchfit)

chart.Histogram(stdret,
     methods = c("add.normal", "add.density"),
     colorset = c("gray", "red", "blue")
)

Leverage effect:

# Specify the GJR GARCH model
garchspec <- ugarchspec(
     mean.model = list(armaOrder = c(0, 0)),
     variance.model = list(model = "gjrGARCH"),
     distribution.model = "sstd"
)

# Estimate the model and compute volatility
gjrgarchfit <- ugarchfit(data = sp500ret, spec = garchspec)
gjrgarchvol <- sigma(gjrgarchfit)

Other specifications:

# Specify AR(1)-GJR GARCH model
garchspec <- ugarchspec(
     mean.model = list(armaOrder = c(1, 2)),
     variance.model = list(model = "gjrGARCH"),
     distribution.model = "sstd"
)

# Estimate the model
garchfit <- ugarchfit(data = sp500ret, spec = garchspec)

# Print the first four coefficients
coef(garchfit)[c(1:4)]
##            mu           ar1           ma1           ma2 
##  0.0004002116  0.6689868348 -0.7027022563 -0.0109852393

APARCH

garchspec <- ugarchspec(
     mean.model = list(armaOrder = c(0, 1)),
     variance.model = list(model = "apARCH"),
     distribution.model = "std"
)
# Estimate the model
aparchfit <- ugarchfit(data = sp500ret, spec = garchspec, solver = "hybrid")
aparchvol <- sigma(aparchfit)
aparchmean <- fitted(aparchfit)

A remark on the skewed Student distribution:

The conditional density to use for the innovations. Valid choices are “norm” for the normal distibution, “snorm” for the skew-normal distribution, “std” for the student-t, “sstd” for the skew-student, “ged” for the generalized error distribution, “sged” for the skew-generalized error distribution, “nig” for the normal inverse gaussian distribution, “ghyp” for the Generalized Hyperbolic, and “jsu” for Johnson’s SU distribution.

# Compare volatility
plotvol <- plot(abs(sp500ret), col = "grey")
plotvol <- addSeries(gjrgarchvol, col = "red", on = 1)
plotvol <- addSeries(garchvol, col = "blue", on = 1)
plotvol <- addSeries(aparchvol, col = "green3", on = 1)
plotvol

GARCH-in-Mean specification and estimation

gim_garchspec <- ugarchspec(mean.model = list(armaOrder = c(0, 0), archm = TRUE, archpow = 2), variance.model = list(model = "gjrGARCH"), distribution.model = "sstd")
gim_garchfit <- ugarchfit(data = sp500ret, spec = gim_garchspec)
gim_garchfit
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : gjrGARCH(1,1)
## Mean Model   : ARFIMA(0,0,0)
## Distribution : sstd 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.000217    0.000073   2.9707 0.002972
## archm   2.038474    0.544626   3.7429 0.000182
## omega   0.000001    0.000002   0.5719 0.567390
## alpha1  0.000740    0.004539   0.1631 0.870441
## beta1   0.913129    0.029509  30.9445 0.000000
## gamma1  0.148437    0.058562   2.5347 0.011255
## skew    0.922223    0.012843  71.8092 0.000000
## shape   7.018567    0.902143   7.7799 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.000217    0.001314 0.165245  0.86875
## archm   2.038474   23.768012 0.085765  0.93165
## omega   0.000001    0.000045 0.027598  0.97798
## alpha1  0.000740    0.023691 0.031248  0.97507
## beta1   0.913129    0.620573 1.471430  0.14117
## gamma1  0.148437    1.243053 0.119413  0.90495
## skew    0.922223    0.145021 6.359238  0.00000
## shape   7.018567   22.134566 0.317086  0.75118
## 
## LogLikelihood : 24374.79 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.6694
## Bayes        -6.6619
## Shibata      -6.6694
## Hannan-Quinn -6.6668
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.3488  0.5548
## Lag[2*(p+q)+(p+q)-1][2]    0.4322  0.7263
## Lag[4*(p+q)+(p+q)-1][5]    4.6044  0.1877
## d.o.f=0
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                      3.022 0.08214
## Lag[2*(p+q)+(p+q)-1][5]     4.002 0.25384
## Lag[4*(p+q)+(p+q)-1][9]     5.162 0.40594
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]     1.359 0.500 2.000  0.2437
## ARCH Lag[5]     1.546 1.440 1.667  0.5803
## ARCH Lag[7]     2.579 2.315 1.543  0.5963
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  397.6714
## Individual Statistics:              
## mu      0.2736
## archm   0.1413
## omega  59.5030
## alpha1  1.2735
## beta1   1.3664
## gamma1  0.5872
## skew    0.8253
## shape   0.9278
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.89 2.11 2.59
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value      prob sig
## Sign Bias            3.234 1.228e-03 ***
## Negative Sign Bias   2.795 5.198e-03 ***
## Positive Sign Bias   2.180 2.929e-02  **
## Joint Effect        28.717 2.568e-06 ***
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     52.19    6.203e-05
## 2    30     54.06    3.184e-03
## 3    40     66.23    4.188e-03
## 4    50     69.03    3.113e-02
## 
## 
## Elapsed time : 0.882395
# Predicted mean returns and volatility of GARCH-in-mean
gim_mean <- fitted(gim_garchfit)
gim_vol <- sigma(gim_garchfit)

The comparison of different conditional means and volatilites

# Correlation between predicted return using one of the previous models and GARCH-in-mean models
cor(aparchmean, gim_mean)
##            [,1]
## [1,] 0.09720179
cor(garch_mean, gim_mean)
##            [,1]
## [1,] -0.1126776
# Correlation between predicted volatilities across models
cor(merge(garchvol, gjrgarchvol, aparchvol, gim_vol))
##              garchvol gjrgarchvol aparchvol   gim_vol
## garchvol    1.0000000   0.9617339 0.9773612 0.9623802
## gjrgarchvol 0.9617339   1.0000000 0.9771208 0.9995686
## aparchvol   0.9773612   0.9771208 1.0000000 0.9773793
## gim_vol     0.9623802   0.9995686 0.9773793 1.0000000

Avoid unneccessary complexity:

# Print the flexible GARCH parameters # flexgarchfit
print(coef(gim_garchfit))
##           mu        archm        omega       alpha1        beta1       gamma1 
## 2.170667e-04 2.038474e+00 1.252247e-06 7.403016e-04 9.131289e-01 1.484370e-01 
##         skew        shape 
## 9.222228e-01 7.018567e+00
# Restrict the flexible GARCH model by impose a fixed ar1 and skew parameter
rflexgarchspec <- gim_garchspec
setfixed(rflexgarchspec) <- list(archm = 1, skew = 1)

# Estimate the restricted GARCH model #EURUSDret
rflexgarchfit <- ugarchfit(data = sp500ret, spec = rflexgarchspec)
rflexgarchfit
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : gjrGARCH(1,1)
## Mean Model   : ARFIMA(0,0,0)
## Distribution : sstd 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.000419    0.000081  5.17139  0.00000
## archm   1.000000          NA       NA       NA
## omega   0.000001    0.000001  1.48500  0.13754
## alpha1  0.001013    0.003950  0.25648  0.79758
## beta1   0.916079    0.010846 84.46270  0.00000
## gamma1  0.144571    0.014695  9.83783  0.00000
## skew    1.000000          NA       NA       NA
## shape   6.719495    0.341236 19.69165  0.00000
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.000419    0.000111 3.761047 0.000169
## archm   1.000000          NA       NA       NA
## omega   0.000001    0.000006 0.190751 0.848721
## alpha1  0.001013    0.024118 0.042008 0.966492
## beta1   0.916079    0.093999 9.745584 0.000000
## gamma1  0.144571    0.149394 0.967716 0.333186
## skew    1.000000          NA       NA       NA
## shape   6.719495    3.063479 2.193419 0.028277
## 
## LogLikelihood : 24360.86 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.6662
## Bayes        -6.6605
## Shibata      -6.6662
## Hannan-Quinn -6.6642
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.5871  0.4435
## Lag[2*(p+q)+(p+q)-1][2]    0.7641  0.5815
## Lag[4*(p+q)+(p+q)-1][5]    5.7898  0.1008
## d.o.f=0
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                      2.680  0.1016
## Lag[2*(p+q)+(p+q)-1][5]     3.700  0.2939
## Lag[4*(p+q)+(p+q)-1][9]     4.884  0.4466
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]     1.334 0.500 2.000  0.2480
## ARCH Lag[5]     1.557 1.440 1.667  0.5774
## ARCH Lag[7]     2.608 2.315 1.543  0.5905
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  430.727
## Individual Statistics:              
## mu      0.2277
## omega  78.1497
## alpha1  1.1330
## beta1   1.2566
## gamma1  0.6079
## shape   0.9506
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.49 1.68 2.12
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value      prob sig
## Sign Bias            3.087 2.029e-03 ***
## Negative Sign Bias   2.544 1.099e-02  **
## Positive Sign Bias   2.168 3.019e-02  **
## Joint Effect        26.943 6.050e-06 ***
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     64.33    7.829e-07
## 2    30     73.14    1.113e-05
## 3    40     88.10    1.164e-05
## 4    50    101.40    1.612e-05
## 
## 
## Elapsed time : 0.7889161
# Compare the volatility of the unrestricted and restriced GARCH models
plotvol <- plot(abs(sp500ret), col = "grey")
plotvol <- addSeries(sigma(gim_garchfit), col = "black", lwd = 4, on = 1)
plotvol <- addSeries(sigma(rflexgarchfit), col = "red", on = 1)
plotvol

Another solution: setbounds()

# Define bflexgarchspec as the bound constrained version
bflexgarchspec <- gim_garchspec
setbounds(bflexgarchspec) <- list(alpha1 = c(0.05, 0.2), beta1 = c(0.7, 0.95))

# Estimate the bound constrained model
bflexgarchfit <- ugarchfit(data = sp500ret, spec = bflexgarchspec)
bflexgarchfit
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : gjrGARCH(1,1)
## Mean Model   : ARFIMA(0,0,0)
## Distribution : sstd 
## 
## Convergence Problem:
## Solver Message:
# Inspect coefficients
coef(bflexgarchfit)
## NULL
# Compare forecasts for the next ten days:
forecast_gim_garch <- ugarchforecast(gim_garchfit, n.ahead = 10)
series_forecast_gim <- forecast_gim_garch@forecast$seriesFor
sigma_forecast_gim <- forecast_gim_garch@forecast$sigmaFor

forecast_garch <- ugarchforecast(garchfit, n.ahead = 10)
sigma_forecast_garch <- forecast_garch@forecast$sigmaFor

cbind(sigma_forecast_gim, sigma_forecast_garch)
##       2017-12-29  2017-12-29
## T+1  0.004782920 0.004617103
## T+2  0.004878952 0.004697538
## T+3  0.004971807 0.004775718
## T+4  0.005061678 0.004851765
## T+5  0.005148739 0.004925786
## T+6  0.005233148 0.004997880
## T+7  0.005315048 0.005068141
## T+8  0.005394570 0.005136651
## T+9  0.005471833 0.005203489
## T+10 0.005546948 0.005268728

#Variance targeting Financial return volatility clusters through time: periods of above average volatility are followed by period of below average volatility. The long run prediction is that:

when volatility is high, it will decrease and revert to its long run average.
when volatility is low, it will increase and revert to its long run average.

In the estimation of GARCH models we can exploit this mean reversion behavior of volatility by means of volatility targeting. We then estimate the GARCH parameters in such a way that the long run volatility implied by the GARCH model equals the sample standard deviation.

# Complete the specification to do variance targeting
garchspec <- ugarchspec(
     mean.model = list(armaOrder = c(0, 0)),
     variance.model = list(
          model = "sGARCH",
          variance.targeting = TRUE
     ),
     distribution.model = "std"
)

# Estimate the model
garchfit <- ugarchfit(data = sp500ret, spec = garchspec)

# Print the GARCH model implied long run volatility
sqrt(uncvariance(garchfit))
## [1] 0.01099808
# Verify that it equals the standard deviation (after rounding)
all.equal(sqrt(uncvariance(garchfit)), sd(sp500ret), tol = 1e-4)
## [1] "Mean relative difference: 0.0003076497"
# Specify model with AR(1) dynamics, GJR GARCH and skewed student t
flexgarchspec <- ugarchspec(
     mean.model = list(armaOrder = c(1, 0)),
     variance.model = list(model = "gjrGARCH"),
     distribution.model = "sstd"
)

# Estimate the model
flexgarchfit <- ugarchfit(data = sp500ret, spec = flexgarchspec)
# Complete and study the statistical significance of the estimated parameters
round(flexgarchfit@fit$matcoef, 6)
##         Estimate  Std. Error    t value Pr(>|t|)
## mu      0.000340    0.000083   4.114964 0.000039
## ar1    -0.029168    0.011651  -2.503390 0.012301
## omega   0.000001    0.000000   2.433509 0.014953
## alpha1  0.001521    0.003511   0.433308 0.664791
## beta1   0.918116    0.007553 121.549276 0.000000
## gamma1  0.143252    0.003340  42.884298 0.000000
## skew    0.917839    0.014330  64.048800 0.000000
## shape   6.932397    0.499364  13.882460 0.000000
# str(flexgarchfit)

A comment for the skew parameter: The skewed student t skewness parameter is statistically significantly different from zero. However 0 is not an interesting value to compare with for the skewness parameter. We should compare it with 1, since the the distribution is symmetric if the skewed student t skewness parameter equals one. Here the estimated skewness parameter is 0.917 and its standard error equals 0.0143. The estimate is thus more than two standard errors away from one and thus we can conclude that the skewness parameter is different from one and that the distribution is thus not symmetric.

# Specify model with constant mean, standard GARCH and student t
tgarchspec <- ugarchspec(
     mean.model = list(armaOrder = c(0, 0)),
     variance.model = list(model = "gjrGARCH", variance.targeting = TRUE),
     distribution.model = "sstd"
)

# Fix the mu parameter at zero
setfixed(tgarchspec) <- list("mu" = 0)

# Estimate the model
tgarchfit <- ugarchfit(data = sp500ret, spec = tgarchspec)

# Verify that the differences in volatility are small
plot(sigma(tgarchfit) - sigma(flexgarchfit))

Comparison of mean errors and variance errors:

# Compute prediction errors
garcherrors <- residuals(garchfit)
gjrerrors <- residuals(tgarchfit)

# Compute MSE for variance prediction of garchfit model
mean((sigma(garchfit)^2 - garcherrors^2)^2)
## [1] 1.293864e-07
# Compute MSE for variance prediction of gjrfit model
mean((sigma(tgarchfit)^2 - gjrerrors^2)^2)
## [1] 1.216576e-07
# Print the number of estimated parameters
print(coef(garchfit))
##           mu       alpha1        beta1        shape        omega 
## 6.732457e-04 7.108829e-02 9.232537e-01 6.377146e+00 6.843804e-07
print(coef(tgarchfit))
##           mu       alpha1        beta1       gamma1         skew        shape 
## 0.000000e+00 1.556731e-06 9.154120e-01 1.533007e-01 9.029445e-01 7.218296e+00 
##        omega 
## 1.311604e-06

or:

# Print the number of estimated parameters
length(coef(garchfit))
## [1] 5
length(coef(tgarchfit))
## [1] 7
# Print likelihood of the two models
likelihood(garchfit)
## [1] 24271.46
likelihood(tgarchfit)
## [1] 24366.18
# Print the information criteria of the two models
infocriteria(garchfit)
##                       
## Akaike       -6.642251
## Bayes        -6.638475
## Shibata      -6.642251
## Hannan-Quinn -6.640952
infocriteria(tgarchfit)
##                       
## Akaike       -6.667901
## Bayes        -6.663182
## Shibata      -6.667902
## Hannan-Quinn -6.666279

Diagnose the standardize absolute returns

# Compute the standardized residuals
stdret <- residuals(tgarchfit, standardize = TRUE)

# Compute their sample mean and standard deviation
mean(stdret)
## [1] 0.03628148
sd(stdret)
## [1] 1.004238
# Compute the standardized returns
stdsp500ret <- residuals(tgarchfit, standardize = TRUE)

# Compute their sample mean and standard deviation
mean(stdsp500ret)
## [1] 0.03628148
sd(stdsp500ret)
## [1] 1.004238
# Correlogram of the absolute (standardized) returns
par(mfrow = c(1, 2))
acf(abs(sp500ret), 22)
acf(abs(stdsp500ret), 22)

# Ljung-Box test
Box.test(abs(stdsp500ret), 22, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  abs(stdsp500ret)
## X-squared = 39.116, df = 22, p-value = 0.01369

Shorter sample:

# Estimate the model on the last 1000 observations
garchspec <- ugarchspec(
     mean.model = list(armaOrder = c(0, 0)),
     variance.model = list(model = "sGARCH"),
     distribution.model = "std"
)
garchfit1000 <- ugarchfit(data = tail(sp500ret, 1000), spec = garchspec)
round(garchfit1000@fit$matcoef, 4)
##         Estimate  Std. Error  t value Pr(>|t|)
## mu        0.0007      0.0002   4.2618    0.000
## omega     0.0000      0.0000   1.2107    0.226
## alpha1    0.2138      0.0478   4.4739    0.000
## beta1     0.7670      0.0453  16.9136    0.000
## shape     4.7618      0.7632   6.2395    0.000
# Compute standardized returns
stdret <- residuals(garchfit1000, standardize = TRUE)

# Do the Ljung-Box test on the absolute standardized returns
Box.test(abs(stdret), 22, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  abs(stdret)
## X-squared = 22.99, df = 22, p-value = 0.4023

In three steps:

# Estimate the GARCH model using all the returns and compute the in-sample estimates of volatility
garchinsample <- ugarchfit(data = sp500ret, spec = garchspec)
garchvolinsample <- sigma(garchinsample)

# Use ugarchroll for rolling estimation of the GARCH model
garchroll <- ugarchroll(garchspec,
     data = sp500ret,
     n.start = 2000, refit.window = "moving", refit.every = 2500
)

# Set preds to the data frame with rolling predictions
preds <- as.data.frame(garchroll)
# Compare in-sample and rolling sample volatility in one plot
garchvolroll <- xts(preds$Sigma, order.by = as.Date(rownames(preds)))
volplot <- plot(garchvolinsample, col = "darkgrey", ylim = c(-0.1, 0.1))
volplot <- addSeries(garchvolroll, on = 1, col = "blue")
plot(volplot, main = "In-sample versus rolling vol forecasts")

# ? not visible

The proof of the pudding is in the eating.

# Inspect the first three rows of the dataframe with out of sample predictions
garchpreds <- as.data.frame(garchroll)
head(garchpreds, 3)
##                      Mu       Sigma Skew    Shape Shape(GIG)      Realized
## 1996-11-29 0.0006421517 0.006061525    0 5.238311          0  0.0026754967
## 1996-12-02 0.0006421517 0.006000423    0 5.238311          0 -0.0006076458
## 1996-12-03 0.0006421517 0.005935000    0 5.238311          0 -0.0109442741
# Compute prediction errors
e <- garchpreds$Realized - garchpreds$Mu
d <- e^2 - garchpreds$Sigma^2

# Compute MSE for the garchroll variance prediction
garchMSE <- mean(d^2)
# Compute MSE for gjrgarchroll
gjrgarchspec <- ugarchspec(
     mean.model = list(armaOrder = c(0, 0)),
     variance.model = list(model = "gjrGARCH"),
     distribution.model = "std"
)

gjrgarchroll <- ugarchroll(gjrgarchspec,
     data = sp500ret,
     n.start = 2000, refit.window = "moving", refit.every = 2500
)

gjrgarchpreds <- as.data.frame(gjrgarchroll)
e <- gjrgarchpreds$Realized - gjrgarchpreds$Mu
d <- e^2 - gjrgarchpreds$Sigma^2
gjrgarchMSE <- mean(d^2)

(gjrgarchMSE - garchMSE) / garchMSE
## [1] -0.04899196
# When doing rolling sample GARCH model predictions, the comparison is fair. You only use the returns available at the time of prediction. As such you are not exposed to look-ahead bias nor to overfitting, since the data used for the prediction is different from the actual return used to evaluate this prediction. The best predicting model wins. Depending on the application, the winner can be a simple model or a complex model.

VaR:

# Extract the dataframe with predictions from the rolling GARCH estimation
garchpreds <- as.data.frame(garchroll)

# Extract the 5% VaR
garchVaR <- quantile(garchroll, probs = 0.05)

# Extract the volatility from garchpreds
garchvol <- xts(garchpreds$Sigma, order.by = time(garchVaR))
# Analyze the comovement in a time series plot
garchplot <- plot(garchvol, ylim = c(-0.1, 0.1), col = "red")
garchplot <- addSeries(garchVaR, on = 1, col = "blue")
garchplot <- addSeries(sp500ret, on = 1, col = "grey")
plot(garchplot, main = "Daily vol and 5% VaR")

# Take a default specification a with a normal and skewed student t distribution
normgarchspec <- ugarchspec(distribution.model = "norm")
sstdgarchspec <- ugarchspec(distribution.model = "sstd")

# Do rolling estimation
normgarchroll <- ugarchroll(normgarchspec,
     data = sp500ret,
     n.start = 2500, refit.window = "moving", refit.every = 2000
)
sstdgarchroll <- ugarchroll(sstdgarchspec,
     data = sp500ret,
     n.start = 2500, refit.window = "moving", refit.every = 2000
)
# Compute the 5% value at risk
normgarchVaR <- quantile(normgarchroll, probs = 0.05)
sstdgarchVaR <- quantile(sstdgarchroll, probs = 0.05)

# Compute the coverage
actual <- xts(as.data.frame(normgarchroll)$Realized, time(normgarchVaR))
mean(actual < normgarchVaR)
## [1] 0.05658415
mean(actual < sstdgarchVaR)
## [1] 0.05887248
# The coverage is clearly better when using a skewed student t distribution. Modelling correctly the distribution is of utmost important for getting the quantiles (and thus the value-at-risk) right. Note that in the implementation we took 2000 observations to reduce the computation time in this example. In practice you will be able to improve performance by re-estimating at a higher frequency.

Fit the model to a particular period:

# Estimate the model
garchfit <- ugarchfit(data = sp500ret["/2006/12"], spec = garchspec)

# Fix the parameters
progarchspec <- garchspec
setfixed(progarchspec) <- as.list(coef(garchfit))

# Use ugarchfilter to plot the estimated volatility for the complete period
garchfilter <- ugarchfilter(data = sp500ret, spec = progarchspec)
plot(sigma(garchfilter))

Compare the forecasts for different periods

garchspec <- ugarchspec(
     mean.model = list(armaOrder = c(0, 0)),
     variance.model = list(
          model = "sGARCH",
          variance.targeting = FALSE
     ),
     distribution.model = "std"
)

# Compare the 252 days ahead forecasts made at the end of September 2008 and September 2017
progarchspec <- ugarchspec(
     mean.model = list(armaOrder = c(0, 0)),
     variance.model = list(model = "sGARCH"),
     distribution.model = "std"
)

### check how it works now:
garchfit2008 <- ugarchfit(data = sp500ret["/2008-09"], spec = garchspec)
garchfor2008 <- ugarchforecast(fitORspec = garchfit2008, n.ahead = 252)
garchfit2017 <- ugarchfit(data = sp500ret["/2017-09"], spec = garchspec)
garchfor2017 <- ugarchforecast(fitORspec = garchfit2017, n.ahead = 252)

# plot
par(mfrow = c(2, 1), mar = c(3, 2, 3, 2))
plot(sigma(garchfor2008), main = "/2008-09", type = "l")
plot(sigma(garchfor2017), main = "/2017-09", type = "l")

# we apply here an earlier specification with GARCH-in-Mean gim_
gim_garchfit <- ugarchfit(data = sp500ret["/2008-09"], spec = gim_garchspec)
garchforecast2008 <- ugarchforecast(data = sp500ret["/2008-09"], fitORspec = gim_garchfit, n.ahead = 252)
gim_garchfit <- ugarchfit(data = sp500ret["/2017-09"], spec = gim_garchspec)
garchforecast2017 <- ugarchforecast(data = sp500ret["/2017-09"], fitORspec = gim_garchfit, n.ahead = 252)

# fitORspec = garchspec
par(mfrow = c(2, 1), mar = c(3, 2, 3, 2))
plot(sigma(garchforecast2008), main = "/2008-09", type = "l")
plot(sigma(garchforecast2017), main = "/2017-09", type = "l")

Simulations

# A code to simulate 5 time series of 10 years of daily returns
spec <- ugarchspec(
     variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
     mean.model = list(armaOrder = c(0, 0), include.mean = TRUE),
     distribution.model = "sstd",
     fixed.pars = list(mu = 0.001, omega = 0.00001, alpha1 = 0.05, beta1 = 0.90, shape = 8, skew = 2)
)

# simulate the path
# path.sgarch = ugarchpath(spec, n.sim=3000, n.start=1, m.sim=1)

simgarch <- ugarchpath(spec = spec, m.sim = 5, n.sim = 10 * 252, rseed = 210)

# Plot the simulated returns and volatility of the four series
simret <- fitted(simgarch)
plot.zoo(simret)

plot.zoo(sigma(simgarch))

# Compute the corresponding simulated prices and plot them
simprices <- exp(apply(simret, 2, "cumsum"))
matplot(simprices, type = "l", lwd = 3)

# Estimate using default starting values
garchfit <- ugarchfit(data = sp500ret, spec = garchspec)

# Print the estimated parameters and the likelihood
coef(garchfit)
##           mu        omega       alpha1        beta1        shape 
## 6.751676e-04 6.167870e-07 7.470453e-02 9.230754e-01 6.087123e+00
likelihood(garchfit)
## [1] 24272.28
# Set other starting values and re-estimate
setstart(garchspec) <- list(alpha1 = 0.05, beta1 = 0.9, shape = 6)
garchfit <- ugarchfit(data = sp500ret, spec = garchspec)

# Print the estimated parameters and the likelihood
coef(garchfit)
##           mu        omega       alpha1        beta1        shape 
## 6.706921e-04 6.143096e-07 7.436183e-02 9.234476e-01 6.110132e+00
likelihood(garchfit)
## [1] 24272.27

#Covariance - with the dataset

# Compute the standardized returns, together with their correlation
stdret_aparch <- residuals(aparchfit, standardize = TRUE)

gim_garchfit <- ugarchfit(data = sp500ret, spec = gim_garchspec)
stdret_gim <- residuals(gim_garchfit, standardize = TRUE)

garchfit_gim_cor <- as.numeric(cor(stdret_aparch, stdret_gim))
print(garchfit_gim_cor)
## [1] 0.9912919

Other series - home:

-download data for MSFT and WMT -fit models -show correlations in the moving window

```{r}

# Compute the covariance and variance of the US and EU returns

msftwmtcov <- sigma(msftgarchfit) * sigma(wmtgarchfit) * msftwmtcor

msftvar <- sigma(msftgarchfit)^2

wmtvar <- sigma(wmtgarchfit)^2

# Compute the minimum variance weight of the MSFT in the MSFT-WMT portfolio

portweight <- (wmtvar - msftwmtcov) / (msftvar + wmtvar - 2 * msftwmtcov)

plot(portweight)

```

```{r}

# Compute standardized returns

stdmsftret <- residuals(msftgarchfit, standardize = TRUE)

stdwmtret <- residuals(wmtgarchfit, standardize = TRUE)

cor(stdmsftret, stdwmtret)

# Load the package PerformanceAnalytics

library(PerformanceAnalytics)

# Plot the 3-month rolling correlation

chart.RollingCorrelation(stdmsftret, stdwmtret, width = 66, main = “3-month rolling correlation between MSFT and WMT daily returns”)

```